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The Gamma distribution is well-known and widely used in many signal processing and communica- 
tions applications. In this letter, a simple and extremely efficient accept/reject algorithm is introduced for 
the generation of independent random variables from a Gamma distribution with any shape parameter 
a > 1. The proposed method uses another Gamma distribution with integer ap < a, from which samples 
can be easily drawn, as proposal function. For this reason, the new technique attains a higher acceptance 
rate (AR) for a > 3 than all the methods currently available in the literature, with AR ^ 1 as a — > oo. 

I. Introduction 
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^ The Gamma probability density function (PDF) is given by Poix) = Cip{x), with Ci = j^^, and 

^H p{x) = x""^ exp (-/3x) , X > 0, (1) 



where a > is the shape parameter, /3 > is the rate parameter and T{a) is the Gamma function |[6l 



^ Chapter 4]. The Gamma distribution is well-known and widely used in different fields, such as Bayesian 

inference [10], signal processing lUTI and digital communications |14|. In particular, in communications 
it has been recently applied in the simulation of fading/shadowing channels using the WeibuU-Gamma 
model H m or the effects of the turbulent atmosphere in free-space optical links with the Gamma-Gamma 
approach, which requires two independent Gamma random variables |[9l [T3]| . 

All of the aforementioned applications need the generation of independent Gamma random variables 
(RVs), X, with arbitrary values of a and /3, i.e., X ~ G{a,(3). When a is an integer, an exact sampler 
is available. Indeed, if a = n G N+ the Gamma PDF becomes an Erlang PDF m, and X can be 



generated as the sum oi n = a independent exponential RVs, i.e., X = Yl^=i -^«' where each Ei follows 
an exponential PDF with parameter /?. These exponentials can be easily obtained through the inversion 
method [6^ Chapter 3], allowing us to express X as 

X = -^J2^nm = -^In (fluA , (2) 

where the Ui are uniform RVs, i.e., Ui ~ Z^([0, 1]). For a ^ N^, the problem of generating a Gamma RV 
X is usually divided in two subcases: a < 1 and a > 1. Focusing on the second case, which is the one 
addressed in this paper, we note that an exact sampler does not exist, but several accept/reject methods 
have been introduced (see (H HI [71 for a review of the approaches proposed). Most of these methods 
consider only /3 = 1, since, given X ~ G{a, 1), it can be easily shown that X = \X ~ G(a,/3). 

In this letter, we develop an extremely efficient rejection sampler to draw independent samples from 
a Gamma PDF with a > 1 and any value of /3. The proposed method outperforms all the alternative 
techniques reported in the literature in terms of acceptance rate (i.e., the key performance measure of a 
accept/reject methods) for a > 3. The main idea is using a suitable Gamma PDF with an integer Up < a 
as a proposal density, from which samples can be easily drawn using Eq. Q. Since the proposal is itself 
another Gamma PDF, it provides a very good fit of the target, thus attaining very high acceptance rates 
that tend to 100% for a — )• +oo, i.e., virtually providing exact sampling. 

II. Background: accept/reject algorithm 

Rejection sampling (RS) is a classical technique for generating independent samples from an arbitrary 
target PDF, Po{x) = Cip{x) with x £ V and Ci = [f^p{x)dx]^^, using an alternative simpler proposal 
PDF, TTo{x) = C2Tt{x) with x ^ V and C2 = [[jy7r{x)dx]^^, such that 7r(2;) > p{x), i.e., 7r(x) is a hat 
function w.r.t. p{x). RS works by generating samples from the proposal density, x' ~ vro(x), accepting 
them when u' < p{x')/7r{x'), with u' ~ U{[0, 1]), and rejecting them otherwise. The key performance 
measure for RS is the average acceptance rate (AR), a/j = J^ ^^7ro{x)dx = ^ < 1. The value of an 
depends on how close the proposal is to the target, and determines the efficiency of the approach. Hence, 
the main difficulty when designing an RS algorithm is finding a good hat function, 7r(x) > p{x), such 
that tt{x) and p{x) are as close as possible and drawing samples from tTo{x) = C2vr(x) can be done 
easily and efficiently. 



III. Novel technique 

In this letter, we consider as target density the PDF given in Eq. ([T]) with a > 1 and any /3 > 0. As 
proposal PDF, we suggest using another Gamma density with different parameters, namely, 

iTo{x) oc tt{x) = Kp x"""^ exp (— /Spx) , x >0, (3) 

where Up = [a\ < a, with [qJ denoting the integer part of q G [1, +oo), and the remaining parameters 
(Pp and Kp) adjusted to obtain the same location and value of the maximum for the proposal and the 
target. On the one hand, for a > 1, 

-' - 1 ^Op-i 



(^P = - = /?^^^. (4) 

max 



^max C^ J- 



K, — 



^ x^p ^ exp {—(ipx) 



a — 1 



(a-ap) 



= exp [Up - a) I — — 1 , (5) 

where Xmax = ^^ is the location of the single maximum of the Gamma PDF, obtained solving -^^ = 0. 
On the other hand, when a = 1, the parameters are set as Op = a, /3p = /3 and Kp = 1. Thanks to this 
choice of Op and the parameters derived in Eqs. (|4]) and Q, we can ensure that: (a) we can draw samples 
exactly from 7ro(x) oc 7r(a;) [6]; (b) it{x) > p{x) for all x > 0, as proved in the following section. 

Our algorithm can be summarized in the following three steps: (1) calculate the parameters of the 
proposal PDF, 'Ko{x) oc 'k{x); (2) draw a sample x' from 7ro(a;) using the direct approach described in 
Eq. Q, i.e., generate Up independent uniform RVs, Uj ~ ^([0, 1]) with i = 1, ..., ap, and set 

x' = --^ln[fjnj ; (6) 

(3) accept x' with probability p{x')/ir{x'), discarding it otherwise. Repeat steps (2) and (3) until the 
desired number of samples have been obtained. 

IV. Proof of the RS inequality 
In order to apply the RS technique, we need to ensure that 7r(x) > p{x), i.e., 

Kpx"'"-^ exp (-/3px) > x"~^ exp {-(3x) , Vx > 0. (7) 

For X > 0, Eq. ^ can be rewritten alternatively as 

X'pexp(f]x) >x"~°^ (8) 



where Q = (5 — (3p and x" "" presents a sub-linear growth, since < a — Op < 1. Finally, taking the 
logarithm on both sides of ([8]l, 

ln{Kp) + nx>{a-ap)ln{x). (9) 

Now, since a > ap and (3p is given by Q, we note that 

1 



Q = /5-/3p = /3M--^l >0. (10) 

Hence, the linear function on the left hand side of (|9]l is increasing. Moreover, since a > ap, the 
logarithmic function on the right hand side of (|9]l is increasing and . Consequently, since both functions 
are increasing and concave for x > (i.e., their second derivatives are lower or equal than zero), they can 
have at most two intersection points. Indeed, in the sequel we show that they are tangent at x = Xmax. 
which is the only contact point between both curves for x > 0. In order to prove this, we need to show: 

(1) that both functions are equal at x = Xmax> i-fi-. 

\n{Kp) + ilxmax = (a - ap) In(xmax) = (a - ap) In f j , (11) 

which is fulfilled by construction of the proposal, as Kp and Pp are set to achieve 7r(xmax) = ^(a^max); 

(2) that their first derivatives are equal, i.e.. 



d{\n{Kp) + Ox) 



dx 
and taking the derivatives we obtain 



d{{a — ap) Inx) 



dx 



(12) 



n = ^^^ = ^-^^^^ = p(i-^], (13) 

^ma,x a 1 \ a 



which is the result given by Eq. (10 1 



Consequently, since x grows faster than ln(x), we can guarantee that \-a{Kp) + ilx > (a — ap) ln(x) 
for X > 0, with equality only at x = Xmax- Hence, the RS inequality in Q is satisfied and the proposal is 
indeed a hat function for p{x), i.e., 7r(x) > p{x), with equality only at x = 0, x = Xmax and x — )• +oo. 

V. Results 

Gamma generators in the literature are usually designed and compared for /3 = 1, without loss of 
generality. Hence, in order to obtain a fair comparison we only consider /3 = 1, although our approach 
is valid for any value of /3. We have compared the acceptance rate (AR), a^, of the different algorithms 
described below ||6l|8ll3|: 



Our method (Ml): The AR can be expressed analytically as aju = /S*^"" ") p/"\ , i.e., for /3 = 1, 
«i?i = r( \ ■ Note that, for a — )■ +00, we obtain uri — )■ 1, i.e., our approach provides exact 
sampling asymptotically. 

^/2^ 



Log-logistic method (M2) |5|: The proposal PDF is 7r(x) 
and Ki = 4a("+^)e"". The theoretical AR is am 
Cauchy method (M3) Oj: The proposal is 7r(x) = K; 
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-A(a- l)('"-i)e^~"', where c = 7r + 2arctan(a/A). The AR is 0^3 
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and we have 



3a-0.75 



and Ko — „ , 

«i?,3 — ^ ^ ~ 0.56, for a — ^ +00. 

T-student method (M4) HI: The proposal is tt{x) = K3(l + 0.5(2^^=^)2)-3/2 ^jth r/ 

and i^3 = (a - l)°~ie^"°. It can be shown that aR4 -> ^ ^ 0.72 for a -^ +00. 

Modified Ratio-of-Uniforms (RoU) (MS) |12|: It is a variant of the RoU scheme mm, relocating 

the mode of p{x) at x = 0. The asymptotic AR is a/j4 — ;■ ^^ ss 0.73, for a — )■ +00. The AR, 0^4, 

is almost constant for all values of a > 1, as shown in Figure [T] 




Fig. 1. Acceptance rate (AR) using our method Ml (continuous line), M2 (squares), M3 (circles), M4 (triangles) and M5 
(x-marks), for 1 < a < 30 and /3 = 1. 



Figure [T] shows the ARs of all the techniques described above, obtained empirically after drawing 
A^ = 6 • 10^ independent samples, for different values of a. For 1 < a < 2 the best technique is M4. 



For 2 < a < 3, the best method can be either M2 or the newly proposed approach (Ml), depending 
on the value of a. For a > 3, our technique (Ml) is extremely efficient, always outperforming the rest 
of the methods and providing the best results ever reported in the literature. Furthermore, our technique 
provides exact sampling (i.e., aju — ;■ 1) asymptotically as a — )■ +oo. 

VI. Conclusions 

We have developed a rejection sampling (RS) scheme for generating Gamma random variables, with 
arbitrary values of a > 1 and /3, where the proposal PDF is itself another Gamma density. The proposed 
algorithm is simple and extremely efficient, providing the best acceptance rates ever reported in the 
literature for a > 3. 
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